Combining measurements from breathing rate sensors

ABSTRACT

A method for combining breathing rate measurements from two or more independent measurement channels. Independent measurements of breathing rate, for instance by impedance pneumography and photoplethysmography, can be combined to derive an improved measurement eliminating artefacts on one channel. A model of the process generating the breathing rate, is constructed and is run independently for each channel to generate predictions of the breathing rate. The model may be a Kalman filter. The measured values are compared with the predicted values and the difference is used as an indication of the confidence in the measurement, the higher the difference the lower the confidence. The measurements from the two channels are combined using weights calculated from the respected differences.

This application is the US national phase of international applicationPCT/GB02/05684 filed 13 Dec. 2002, which designated the US.PCT/GB02/05684 claims priority to GB Application No. 0130010.2 filed 14Dec. 2001. The entire contents of these applications are incorporatedherein by reference.

BACKGROUND AND SUMMARY

The present invention relates to a method and apparatus for measuringthe breathing rate of a subject, such as a human or animal, and inparticular to a way of combining measurements from two or more breathingrate sensors in order to provide an improved measurement of breathingrate.

There is no clinically acceptable method for the non-invasivemeasurement of breathing rate. A nasal thermistor provides a simple andinexpensive means of tracking respiration but it is obtrusive and is notdeemed to be acceptable for patients on the general ward or even theCoronary Care Unit (CCU). Electrical impedance plethysmography (alsoknown as impedance pneumography) can be used to provide an indirectmeasure of respiration by measuring the changes in electrical impedanceacross the chest with breathing. In this method a small-amplitude,high-frequency current is injected into the body through a pair ofsurface electrodes and the resulting voltage is demodulated to obtainimpedance measurements. The electrical impedance increases ashigh-resistivity air enters the lungs during inspiration but part of thechange is also due to the movement of the electrodes on the chest wall.When monitored in a clinical environment, the impedance plethysmography(IP) signal is often very noisy and is seriously disrupted by patientmovement or change in posture. As a consequence, it has not beenconsidered reliable enough to provide respiration information forregular use on the ward.

Respiratory information is found in other signals recorded from patientswith non-invasive sensors. For example, both the electrocardiogram (EGG)and photoplethysmogram (PPG) waveforms are modulated by the patient'sbreathing. The PPG signal represents the variation in light absorptionacross a finger or earlobe with every heart beat. This signal ismeasured at two wavelengths (usually in the red and near infra-red partsof the spectrum) in a standard pulse oximeter.

FIG. 2 shows 30-second sections of IP (FIG. 2A) and PPG (FIG. 2B)signals recorded from a patient in CCU. Referring first to the IP signalin FIG. 2A, the main peaks I, each corresponding to a breath, can beidentified. However, changes in electrical impedance with the heart beat(as opposed to respiration) are apparent between t=1607 and t=1608(marked as “A”) and probably conceal a peak caused by a breath.Referring to FIG. 2B, the modulation of the PPG caused by respiration(marked “I”) can be seen, although there are small movement artefacts B,C at t=1617 and t=1629. FIG. 3 shows the IP and PPG signals from thesame patient, for a five-minute period. The effect, marked D, E, ofmovement artefact on the finger probe from which the PPG signal isrecorded becomes very obvious at t=1770 and t=1920.

Clearly the artefacts caused by movement and heartbeat would affect themeasurement of breathing rate from these signals. One might considerremoving these artefacts by some type of filtering and thresholding.

FIG. 4B shows the effect of applying a Finite Impulse Response (FIR)low-pass filter with a pass-band cut-off at 0.33 Hz (corresponding to 20breaths per minute) and a stop-band cut-off of 1 Hz (−50 dB) to the IPwaveform of FIG. 4A. The cardiac-synchronous changes are filtered out(see the output of the filter between t=1605 and t=1610) and therespiratory cycle is clear. There are as many peaks (i.e. breaths) inFIG. 4B as there are in FIG. 4D, which shows the result of tracking thepeaks of the PPG signal shown in FIG. 4C and interpolating (withstraight lines) between each peak. The modulation envelope picked out bythe peak tracking produces another respiratory waveform. The sameinformation is again displayed on a five-minute timescale in FIG. 5. Thebreathing rate can be estimated by calculating the interval between twosuccessive peaks of these waveforms, inverting the result andmultiplying it by a factor of 60 in order to obtain an estimate ofbreathing rate in breaths per minute.

The results of these computations over the five-minute period are shownin FIGS. 5E and 5F respectively. The breathing rate is approximately 18breaths/minute throughout the period, but the occurrence of peaks causedby cardiac as opposed to respiratory changes in the IP signal of FIGS.5A and 5B gives rise to erroneously high breathing rates at t=1720,1730, 1780, 1795, 1830 and 1860. Unfortunately no amount of optimisationof the FIR low-pass filter characteristics will ever completely removethe cardiac-synchronous information which is occasionally veryprominent. For instance, in a hyperventilating patient, the breathingrate may be as high as 40 breaths per minute, which is similar to a slowheart rate. So separating signal from noise on a fixed basis isimpossible.

Similarly, the two major instances D, E of movement artefact in the PPGsignal of FIGS. 5C and 5D at t=1770 and t=1920 cause erroneously lowestimates of breathing rates at D′ and E′ because there is a significantdelay between the last valid peak and the first peak after the movementartefact. The estimate of breathing rate from the PPG signal is alsomore variable because occurrences of even slight movement artefactaffect the tracking of the peaks of the PPG signal.

The present invention provides a method and apparatus for improvingmeasurement of breathing rate by combining two measurements of it in away which allows valid changes in the breathing rate to be distinguishedfrom artefacts. In accordance with the invention two measurements ofbreathing rate made in different ways are combined with weights based onthe amount of “confidence” in the measurement, to give an improvedmeasurement or estimate of the actual breathing rate. The twomeasurements may, for example, be obtained using impedance pneumographyand photoplethysmography, though other signals which include respiratoryinformation (for instance ECG) can also be used instead of either ofthese signals, or in addition to them.

In more detail the invention provides a method of measuring breathingrate of a subject comprising the steps of: predicting the value of eachof two independent measurements of the breathing rate, making twoindependent measurements of the breathing rate to produce two measuredvalues, calculating the respective differences between the predictedvalues and the measured values, and combining the two measured valueswith weights determined by said differences.

The steps of prediction, measurement, calculation and combination may berepeated continuously. The predicted value for each of the independentmeasurements may be based on the preceding predicted value and thedifference between the preceding predicted value and the precedingmeasurement.

The predicted value for each of the measurements may be calculated usinga linear or non-linear predictive model, and the model may be adaptive,adapting in dependence upon the amount of process noise in themeasurements.

In the combining of the two measured values, the weight of each valuemay vary inversely with the square of the difference between thepredicted value and the measurement. In one example the two measuredvalues may be combined according to the formula:

${BR} = {{{BR}_{1}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}} + {{BR}_{2}\frac{\sigma_{1}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}}}$where BR₁ and BR₂ are the two measured values of breathing rate and σ₁and σ₂ are the differences between the two measured values and theirrespective predicted values.

The predicted values for the respective measurements may be based onrespective models of the system, and the models may include estimatesfor process noise and sensor noise. The respective models may bemutually independent and may include the same estimates for processnoise and sensor noise. The models may be Kalman filters.

Measurements for which the differences between both measured values andtheir predicted values exceed a predetermined threshold may bediscarded. Further, artefacts, for instance caused by movement orheartbeat in the measurements may be identified based on the values ofthe differences between the measured values and their predicted values.This identification may be used to discard sections of the signal.

Thus with the present invention a prediction is made for each breathingrate measurement and the actual measurement is compared with itsprediction. The difference is computed, which is termed the“innovation”, and this innovation is used to calculate a weight whichwill be given to that measurement when it is combined with the othermeasurement, also weighted according to its innovation. The weights arecalculated so that if the innovation on one measurement channel is high,whereas the innovation on the other measurement channel is low, themeasurement from the low innovation channel is more heavily weighted.This is because a high level of innovation from one channel coincidingwith a low innovation on the other channel is regarded as indicative ofan artefact on the higher innovation channel.

BREIF DESCRIPTION OF THE DRAWINGS

It will be appreciated that the invention can be embodied using computersoftware, and thus the invention extends to a computer program forcontrolling and executing the method, or parts of it, and to a computerreadable storage medium carrying the program. The invention also extendsto corresponding apparatus for carrying out the method.

The invention will be further described by way of non-limitative examplewith reference to the accompanying drawings in which:—

FIG. 1 schematically illustrates apparatus in accordance with anembodiment of the invention connected to a patient;

FIGS. 2A and 2B show 30-seconds of IP and PPG signal respectivelyrecorded from a patient;

FIGS. 3A and 3B show five-minutes worth of IP and PPG data respectivelyfrom the same patient as FIGS. 2A and B;

FIGS. 4A to D show the IP and PPG data of FIGS. 2A and B and processedversions of that data;

FIGS. 5A to F show IP and PPG data corresponding to FIG. 3, processedversions of that data and breathing rate estimates based on that data;and

FIGS. 6A to G show the same data as FIGS. 5A to F and in addition inFIG. 6G a breathing rate estimate obtained by an embodiment of thepresent invention.

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS

An embodiment of the invention will now be described in which theinvention is applied in the medical field for the measurement ofbreathing rate using the two signals described above, namely theimpedance pneumography and photoplethysmography signals. The results ofapplying this to the data of FIGS. 2 to 5 will be illustrated.

FIG. 1 schematically illustrates a breathing rate measurement apparatusin accordance with an embodiment of the invention. A PPG signal isobtained, as illustrated from a finger of a patient, though the signalcan also be obtained from the earlobe, using a conventional PPG sensor 1which is driven by and supplies its output to PPG apparatus 3. An IPsignal is obtained by conventional IP equipment 7 using two electrodes 5(though four electrodes may be used, with a pair for current injectionand a pair for voltage measurement). The IP apparatus 7 and PPGapparatus 3 supply their signals to a processor 10 which processes thesignals in the manner to be described below, and outputs a display ofbreathing rate on display 12. The display may be in graphical or numericform, and may be combined with a display of others of the patient'svital signs. It will be appreciated that in this embodiment the IPapparatus 7 and PPG apparatus 3 produce an output to the processor 10 ofraw signals of the types illustrated in FIGS. 6A and C. The processor 10processes these signals to produce the signals of FIGS. 6B and D to G.However, the functions of the parts of the apparatus shown in FIG. 1 maybe combined in different ways.

The processor 10 in essence runs a Kalman filter on each breathing ratesignal, and then fuses the filtered signals to produce a fused estimateof the breathing rate.

The Kalman filter is a generic framework for analysing linear dynamicalsystems. Using the process model, the next state x_(t+1) is computedfrom the current state x_(t) using the transition matrix A, assumingfirst-order (Markov) dynamics. The observation model relates theobservations y₁ to the state x₁ of the system via the observation matrixC. The process and observation noise are assumed to be independent andto have zero-mean, Gaussian probability distributions. In the normalphysiological state, the breathing rate can be modelled with a scalarKalman filter which assumes that the time to the next breath is the sameas the interval since the previous breath plus some process noisecharacterising normal or breath-to-breath variability. Thus:x _(t+l) =Ax _(t) +w   (Process model)y _(t) =Cx _(t) +v   (Observation model)

In this embodiment it is assumed that C=1, i.e. the breathing rate isboth the measurement and the state describing the process, with v as thesensor noise model. It is also assumed that the next breathing rate onlyvaries by a small amount with respect to the current one, i.e. A=1 and wis the variance associated with breath-by-breath variability.

A scalar Kalman filter is run separately on each breathing rate signal,using the same process and measurement noise models for both channels.For each channel, the next state x_(pred) is predicted using the Kalmanfilter equation (in which the next state is equal to the previous stateplus the Kalman gain times the innovation). In this embodiment, the nextmeasurement y_(pred) is the same as the next state x_(pred) (since C=1)and, after the next measurement of breathing rate, y_(t+1), theinnovation ε_(t+1) is computed, where:y _(t+1) =y _(pred)+ε_(t+1)ε_(t+1) the difference between the actual value and the predicted value,should normally be a zero-mean white noise sequence. The square of theinnovation, or variance, σ_(t+1) ²=ε_(t+1) ², is the inverse of the“confidence” which is associated with the prediction. A robust estimateof the breathing rate is now obtained by mixing the two breathing rateestimates in inverse proportion to the variance associated with eachone:

${BR} = {{\frac{\sigma_{2}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}{BR1}} + {\frac{\sigma_{1}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}{BR2}}}$where BR1=y_(t+1) and σ₁ ²=σ_(t+1) ² for channel 1 (filtered IP), withequivalent expressions for BR2 and channel 2 (PPG peaks). The innovationvariance and the Kalman gain are also calculated to update the stateestimate in readiness for the next cycle. An example of animplementation of this model in MATLAB is given in Appendix 1. Thatexample is general and will work for vector quantities though in thisembodiment the quantities are scalar. It can be seen from Appendix 1that the predicted breathing rate for each new measurement cycle (xnew)is equal to the previously predicted value (xpred) plus the Kalman gainK times the innovation e. The Kalman gain K is derived from thepredicted variance Vpred and the measurement variance R. The predictedvariance is derived from the previous predicted variance and the processnoise Q. To start the process off it is initialised using an initialvalue of the breathing rate as 15 and an initial value of the statevariance of 40. The process noise in the Q in this embodiment is set to10 and the measurement noise variance R is set to 100.

It will be clear from the implementation that, as normal with a Kalmanfilter, the variance and Kalman gain are not dependent on themeasurement values. The measurement values are only used in the newprediction of breathing rate via the innovation e. Thus it will be notedthat for the constant values of Q and R used in this example the Kalmangain K tends to about 0.2 and the state co-variance V tends to about 27.However, K can be made adaptive by modifying the values for the varianceconstants. Q and R, preferably the process variance Q, according to thetype of process being encountered.

The advantages of this method can be appreciated from four possiblecontexts:

-   1. normal breathing: low innovations on both channels; both    measurements BR1 and BR2 are weighted equally.-   2. valid change seen on both channels: the subject begins to breathe    more quickly or more slowly due to a physiological change; although    the innovation is high, it is high for both channels, and so both    (valid) measurements are again weighted equally.-   3. artefact on one channel: high innovation on one channel only; the    information from that channel is ignored because it is given a low    weighting (high variance).-   4. artefact on both channels: the information is corrupted on both    channels. Prolonged movement artefact, however, is characterised by    high values of innovation on both channels for a sustained period of    time and this can be the basis for discarding sections of data    corrupted by movement artefact.

Note that a fused estimate is computed every time a new measurement ofBR1 or BR2 becomes available, i.e. twice during the respiration cycle(once for each sensor).

It should also be noted that in the event of a problem with an electrodesuch as it becoming detached or falling off, or of signal loss throughsome other cause, this would result in a high innovation and thus a lowweighting for that channel.

FIG. 6G illustrates the result of combining the breathing rate signalsof FIGS. 6E and F in this way.

FIG. 6G illustrates the result of applying the invention to thebreathing rate signals of FIGS. 6E and F obtained respectively from theimpedence pneumography and PPG signals of FIGS. 6A to D. Thus theoriginal IP signal of FIG. 6A is filtered using a finite impulseresponse low-pass filter to produce the processed signal of FIG. 6B. ThePPG signal of FIG. 6C is subjected to peak tracking as described aboveto produce a respiratory waveform of FIG. 6D. Each of these is used toestimate a breathing rate, by calculating the interval between twosuccessive peaks, inverting the result and multiplying it by a factor of60. These two breathing rate estimates are shown in FIGS. 6E and F. Thetwo breathing rate estimates BR1 and BR2 are subjected to the Kalmanfiltering and fusing the process described above and this gives thefused breathing rate estimate of FIG. 6G. It can be seen that apparentchanges in breathing rate which appear only on one channel (for instancemarked F on FIG. 6E or the changes D′ and E′ on FIG. 6F) are removedfrom the fused breathing rate estimate. However, a change in breathingrate on both channels, such as marked G in FIGS. 6E and 6F, appears asG′ on the fused estimate of FIG. 6G.

APPENDIX 1 load -ascii ip_resprate; data_file = ip_resprate; time =data_file(:,1); br = data_file(:,2); start = 1 stop =size(data_file),1); fprintf(‘number of breaths detected = %d \n’, stop);br_limit = 40; % X(t+1) = A X(t)+noise(Q) - process model with Q asvariance of noise w % Y(t) = C X(t)+noise(R) - measurement model with Ras variance of noise ν ss = 1; % state size - sets to one dimensional,ie scalar though routine works for vectors os = 1; % observation size -sets to one dimensional, ie scalar A = [1]; % assume x(t+1)=x(t) C =[1]; % assume y=x Q = 10.0*eye(ss); % process noise - eye is theidentity matrix in MATLAB -here just unity R = 100.0*eye(os); %measurement noise variance initx = [15]; % initial state value(BR of 15bpm) initV = 40*eye(ss); % initial state variance xnew = initx; % -initialisation Vnew = initV; no_update = 0; n = 1; for i = start:stop% - start of cycle x = xnew; % update from previous cycle V = Vnew; %update from previous cycle xpred = A*x; % prediction of state Vpred =A*V*A′ + Q; % prediction of state covariance, A′ is transpose of Aypred(i) = C*xpred; % prediction of measurement y(i) = br(i); % “makemeasurement” t(i) = time(i); % if y(i) == br_limit no_update = 1;%artefact % else no_update = 0; % end e = y(i) − ypred(i); % calculateinnovation innov(i) = e; % for plotting sigma2(i) = e*e; % variance forsaving S = C*Vpred*C′ + R; % innovation covariance Sinv = inv(S); %invert S to compute Kalman gain K = Vpred*C′ *Sinv; % Kalman gain matrixif no_update == 0 % only update if no artefact xnew = xpred + K*e; %update state by the innovation controlled by the Kalman gain Vnew =(eye(ss) − K*C)*Vpred;  % update state covariance end if i − (n*10) == 0fprintf( {grave over ( )}\n %d0 breaths processed {grave over ( )}, n);n = n + 1; end end

1. A method of measuring breathing rate of a subject comprising thesteps of: predicting a value of each of two independent measurements ofthe breathing rate for two independent measurement channels, making twoindependent measurements of the breathing rate via said two independentmeasurement channels to produce two measured values of the breathingrate, calculating respective differences between the predicted valuesand the measured values for each of said two independent measurementchannels, and combining the two measured values with weights determinedby said differences.
 2. A method according to claim 1 in which the stepsof prediction, measurement, calculation and combination are repeatedcontinuously, the predicted value for each of the two independentmeasurements being based on the preceding predicted value and thedifference between the preceding predicted value and the precedingmeasurement.
 3. A method according to claim 2 in which the predictedvalue for each of the two independent measurements is calculated byusing a linear predictive model.
 4. A method according to claim 2 inwhich the predicted value for each of the two independent measurementsis calculated by using a non-linear predictive model.
 5. A methodaccording to claim 2 in which the predicted value for each of the twoindependent measurements is calculated by using an adaptive predictivemodel, and wherein the model adapts in dependence upon the amount ofprocess noise in the measurements.
 6. A method according to claim 1 inwhich in the step of combining the two measured values the weight ofeach value varies inversely with the square of the difference betweenthe predicted value and the measurement.
 7. A method according to claim6 in which the two measured values are combined according to theformula:${BR} = {{{BR}_{1}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}} + {{BR}_{2}\frac{\sigma_{1}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}}}$where BR₁ and BR₂ are the two measured values, and σ₁ and σ₂ are thedifferences between the two measured values and their respectivepredicted values.
 8. A method according to claim 1 in which thepredicted values for the respective measurements are based on respectivemodels of the system.
 9. A method according to claim 8 in which themodels of the system include estimates for process noise and sensornoise.
 10. A method according to claim 8 in which the respective modelsof the system are mutually independent.
 11. A method according to claim10 in which the respective models of the system include the sameestimates for process noise and sensor noise.
 12. A method according toclaim 8, in which the respective models of the system are Kalmanfilters.
 13. A method according to claim 1 further comprising the stepof discarding series of measurements for which the differences betweenboth measured values and their predicted values exceed a predeterminedthreshold for a predetermined period of time.
 14. A method according toclaim 1 in which the two independent measurements are made by impedancepneumography and photoplethysmography.
 15. A method according to claim 1in which there are more than two measurements.
 16. A method according toclaim 1 further comprising the step of identifying artefacts in themeasurements based on the values of the differences between bothmeasured values and their predicted values.
 17. A computer programencoded on a computer-readable medium and comprising program code which,when executed, performs the method of claim
 1. 18. Apparatus formeasuring breathing rate of a subject, the apparatus comprising: twoindependent measurement channels for receiving two independentmeasurements of a breathing rate as two measured values of the breathingrate; and a processor configured to predict values of each of the twoindependent measurements of the breathing rate for the two independentmeasurement channels, to calculate respective differences between thepredicted values and the measured values for each of said twoindependent measurement channels, and to combine the two measured valueswith weights determined by said differences.
 19. The apparatus accordingto claim 18, further comprising: a display for displaying the breathingrate.
 20. Apparatus for measuring breathing rate of a subject, theapparatus comprising: first and second independent measurement channelsfor providing respective independent measurements of a breathing rate asfirst and second measured values of the breathing rate; and a processorconfigured to predict values of each of the first and second independentmeasurements of the breathing rate for the respective first and secondindependent measurement channels using respective Kalman filters, tocalculate respective differences between the predicted values and themeasured values for each of the first and second independent measurementchannels, and to combine the first and second measured values withweights determined by the differences.